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ABSTRACT 


Felix Bloch’s 1928 article made a prediction concerning the dynamical behavior 
of electrons in a solid, subject to a uniform, static electric field. This aspect of his WOIK, 
as later clarified by Zener, showed that electrons accelerated by an electric field in a 
periodic potential, under the right conditions, would oscillate. A theoretical debate as to 
the existence of this phenomenon has been ongoing since Bloch’s proposal. One of the 
most controversial consequences of this prediction is that an electron undergoing Bloch 
oscillations would radiate. The controversy on the theoretical analysis was due to the 
great difficulty in systematically and reliably treating interband transitions by analytical 
methods based on the time-dependent Schrodinger equation for independent electrons. In 
this thesis, we numerically solve the time-dependent Schrodinger equation to show that 
electrons accelerated by in an electric field in periodic structures do undergo Bloch 
oscillations and other dynamic behavior. By accurately modeling this phenomenon, we 
hope to gain a better understanding of it in hopes of using it in future applications as a 


stable source of Terahertz (THz) radiation. 
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L INTRODUCTION 


In 1928, Felix Bloch [Ref. 1] wrote a groundbreaking work discussing the 
quantum mechanical nature of electrons in solids. In this seminal work, Bloch 
demonstrated what is widely known as “Bloch’s theorem,” which establishes the 
mathematical form of the electron wave function of a crystalline solid. We will discuss 
Bloch’s theorem later in this thesis. Bloch’s theorem is the basis for our understanding of 
the all-important energy bands in solids. In addition to establishing Bloch’s theorem, a 
less well known result from Bloch’s 1928 article is a prediction concerning the dynamical 
behavior of electrons in a solid, subject to a uniform, static electric field. This aspect of 
his work, as later expanded by Zener [Ref. 2], showed that electrons accelerated by an 
electric field in a periodic potential, under the nght conditions would oscillate. The 


frequency of oscillation is given by the formula @, = eFa/h , where eis the charge of an 


electron, F is the electric field amplitude, ais the lattice repeat distance, and 1 is 
Planck’s constant h divided by 2x. This phenomenon has come to be known as Bloch 
oscillations. A theoretical debate as to the existence of this phenomenon has been 
ongoing since Bloch’s proposal. One of the most controversial consequences of this 
prediction is that an electron undergoing Bloch oscillations would radiate. The 
controversy on the theoretical analysis was due to the great difficulty in systematically 
and reliably treating interband transitions by analytical methods based on the time- 
dependent Schrodinger equation for independent electrons. 

With the advent and refinement of semiconductor technology, it has been finally 
possible to test the predictions of Bloch’s proposal. In 1970, Esaki and Tsu [Ref. 3] 
proposed using 2 semiconductor superlattice, which would provide the needed periodic 
structure on a large enough scale in order to conduct tests to generate and detect Bloch 
oscillations. According to Esaki and Tsu [Ref. 3], the oscillating electron ina 
semiconductor superlattice would provide a source of Terahertz (THz) radiation. It was 
only in 1992, however, that Bloch oscillations and the associated Terahertz radiation were 


conclusively detected in semiconductor superlattices [Ref. 4]. These recent four wave 
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mixing experiments confirm the presence of photons emitted as a result of the oscillating 
electron, but noted the oscillations died between 1 and 15 oscillations [Ref. 4]. More 
recently, Luban and co-workers[Ref. 5] modeled the effect of an electron in a 
semiconductor superlattice using computer techniques, and showed Bloch oscillations 
exist at Terahertz frequencies, but did not show the decay as seen in real world 
experiments. The advantage of a purely numerical approach to solving the time- 
dependent Schrodinger equation is that one explicitly works with the full Hamiltonian, 
and not some approximation as in analytical solutions, thus including all interband 
transitions in the simulation. 

This paper will cover the theoretical debates and arguments surrounding the issue. 
Also, it will develop the theory behind Bloch oscillations. Lastly, it will cover the 
numerical techniques involved in solving the Time Dependent Schrodinger Equation, and 


present the results of the simulations. 


Il. THEORETICAL WORK 


In this chapter, we summarize and discuss the works of major contributors on the 


subject of Bloch oscillations. 


A. FELIX BLOCH 


In 1928, Felix Bloch wrote a ground breaking paper [Ref. 1] which until recently 
was the subject of great debate. In it, he used quantum mechanics to model the effects of 
electrons in periodic (atomic) structures. In this landmark paper, he proved what is now 
known as Bloch’s theorem, which shows that the mathematical form of the eigenvector of 


the time independent Schrodinger equation must be: 


We (X) = eu, (X) (2.1) 
where 
U, ,(X) = U,, (x +a) (2.2) 


and n labels the energy bands and is known as the band index. Also, the parameter k is 
tne allowed wavevector of the electron and 1s restricted to the Brillouin zone associated 
with the particular lattice structure of the crystal. The associated energy eigenvalue 
E(k) is called the band structure function. 

Bloch showed [Ref. 1] that the wave function of any particle in a periodic 
potential can be represented as a sum of orthogonal functions that have the periodicity of 
the lattice structure. Bloch [Ref. 1], starting with the time dependent Schrodinger 
equation with an external electric field applied: 

2 i2m Ow(x,t) | 


(V(x) — eFxy (x,t) + ———> 0, (2.3) 
nm h ot 





V-w(x,t) - 


and assuming the following form of the wave function: 


2ni 


-—E,t 
BOSE) = Bene nn? ai: (2.4) 


derived the result for the time variation of the amplitude of the coefficients of the Bloch 


functions of the electron wave function: 


2neF oO 5 
—lc_,(0)°. Dud 
h dk | Gnas ( } ( ) 


d 
—|c,,\° = - 


at 





He also showed for a collection of electrons obeying Fermi statistics [Ref. 1], that the 


time variation of the electron distribution function f(k,t) is: 


df  _ 2neF of 


2.6 
at hn OK ea) 


In (2.6), f(k,t) is the probability of finding an electron between states k and k + Ak at 


time ¢. Bloch claims this is the same result obtained by Lorentz [Reif. 6] in an earlier 


paper published in 1916. 


B. CLARENCE ZENER 


Zener [Ref. 2], in 1933, substantially clarified Bloch’s work [Ref. 1] in a paper 
investigating the breakdown of solids by an electric field. The breakdown of a dielectric 


under the influence of an electric field takes place by one of two mechanisms [Ref. 2}: 


1. A process similar to the electrical breakdown of a gas, i.e., avalanche 
breakdown. 


2. A process like the auto ionization of free atoms in an electric field. 


For solids, Zener [Ref . 2] said the second process is more important in a solid than ina 


gas since the mean free path of an electron is much less for a solid than a gas. An 
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electron, therefore, tends to become excited and move about in the lattice instead of being 
freed from the solid completely. He used the Bloch model of an electron to calculate the 
rate at which electrons move from lower to upper energy bands in a solid under the 
influence of a constant electric field. 

Zener [Ref. 2] started with the same initial wavefunction (2.1). He said that the 
condition that w, be finite over the lattice imposes the restriction that k be real as well. 
This only occurs when the electron energy E, lies within the allowed energy bands that 
are a Solution of the time independent Schrodinger equation, with the first energy band 
occuring within the values -x/a s k < x/a[Ref. 2] for a one dimensional model. 
Following Bloch, Zener [Ref. 2] expanded his wave function in terms of wave functions 


of the lowest energy band with no electric field applied: 


a/a 
WEE) =f BD na Mk 2.7) 


Following Bloch’s derivation, Zener [Ref. 2] reproduced equation (2.5) in a more concise 





form: 
ae 2meF cn 5 / 
hee __, Aisling yess 2.8 
=“ le--— ale 2.8) 
where 
lal’ = ok = al ; ; (2.9) 
he 


with G being an arbitrary function. Zener [Ref. 2] also showed that the velocity in 


k space is 2neF/h. Using the limits of k , Zener showed the period of time for real 


space oscillation is: 


T, = (2n/a)/(2neF/h) = — (2.10) 


or the frequency of oscillation is given by: 
0, S==—. (2a ) 


Zener [Ref. 2] said than an electron confined to a single energy band will move opposite 
to the direction of the field until reflected by the lattice, and then it moves in the opposite 
direction until it is stopped by the field, where it starts the same motion over again. 

Zener [Ref. 2] then calculated the probability that an electron would undergo transitions 
between energy bands. Starting with a wave function that is periodic in time, Zener 
showed the wave equation for an electron when an external field F is present is (2.3). 
Zener [Ref. 2] found an approximate solution, assuming the term eF*x varies slightly over 


the lattice distance a: 


Way(x) =e! u, (x) (2.12) 


where Kis a function of (EF - eFx). Using (2.12), and solving for K , Zener [Ref. 2] 


showed the rate of transition (y ) between energy bands is given by: 





eFa x? maE* 
= ——exp| —-— 2A 


where Eis the energy gap between bands. 


C. W. V. HOUSTON 


In 1939, Houston [Ref. 7] showed that the wave vector of an electron accelerated 
by an electric field in a periodic potential grows linearly with time within the bounds of 


the first Brillouin zone[Ref. 7]. 


Houston assumed the form of this wave function as: 


{ 
~$ [Eksnede 


W(X, t) = Ua (XE (2.14) 


where A = eF /h [Ref. 7], which shows the solution of the time dependent Schrodinger 
equation is constructed from the solution of the time independent Schrodinger equation, 
modulated by a time dependent wave vector. Equation (2.14) is a solution of (2.3), with 


the following remainder: 


f 
oat Ex+)24t 


-(F\ 0%, lee ae (2215) 
[ 


The remainder term is zero for the case when no electric field is present (A = Q), and is 
zero when Vu,,,, = 0 for the free electron case [Ref. 7]. In the case where Aor Vu,,,, are 
not zero, but their product is very small, Houston [Ref. 7] claimed (2.14) is a good 
approximate solution to (2.3). In the case where the remainder is not small, Houston 

[ Ref. 7] showed that the electron wave function can be written as an infinite sum of 


orthogonal Bloch eigenfunctions: 


t 
i = 
hf Fk+221 +t 


i¢k+2n1 +22 
~(x,t) = > a,(t)u, 7, (x)e- e (2.16) 


where / are reciprocal lattice vectors associated with the crystal structure. 


~) 


Using the method of variation of constants, Houston solved for the time variation of the 


coefficients a,(t). Having solved this, he showed [Ref. 7] that the probability for the 


electron to transition across a Brillouin zone boundary is: 


P=4n7e*/A+e%) (2.17) 


where a = mV°a? /eFh* . Houston claims [Ref. 7] this result is similar to Zener’s [Ref. 


2). 
D. GREGORY WANNIER 


Wannier’s contribution to the field of solid state physics is immense. One of his 
many contributions was his work on the motion of electrons in solids in the presence of a 
uniform electric field, and the nature of the energy bands formed by the electric field in 
that solid. Ina series of papers from 1960 to 1968 [Refs. 8-11], Wannier decribed the 
effects of electric and magnetic fields on electrons in periodic potentials. Wannier’s [Ref 
.8] most significant contribution was to show, for a single band Hamiltonian, that the 
energy levels of a solid under the effect of an electric field become discretized such that 


they are a function of the lattice repeat distance a: 


x/a 
E, = [W,(k)dk + nFa (2.18) 
0 


where [Ref. 9] EF are the energy levels, the integral term is the mean energy of the band, 


and nFa are the ladder spacing intervals. This series of equally spaced energy levels is 
termed the Stark ladder, or sometimes the Wannier-Stark ladder. The basic idea behind 
the Stark ladder is that the electric field destroys the degeneracy of the atomic energies of 
the solid that form the bands. Thus, in Wannier’s picture, electric field destroys the zero- 
field energy bands and replaces them with a ladder of equally spaced energy levels. The 


existence of Stark ladders in solids has itself been as controversial a topic as Bloch 


oscillations, and Stark ladders were first observed in semiconductor superlattices for a 
single band beginning in the 1980’s [Ref. 12]. Wannier also showed [Ref. 10] that if the 
electron wave function is comprised of Bloch functions, then the energy bands are 
decoupled, which supports his idea of a discrete energy bands. 

Perhaps Wannier’s most widely known contribution to solid state physics are the 
“Wannier functions,” which are in essence a Fourier transform of the Bloch states [Ref. 
13]. Whereas Bloch states are extended states, i.e., the electron can be found over the 
entire crystal, the Wannier functions are localized to within a few lattice sites. It can be 
shown that Wannier functions centered about different lattice sites are orthogonal [Ref. 
13]. Also, it was proved that Wannier functions at the same lattice site, but representing 
different energy bands are orthogonal [Ref. 13]. The Wannier functions are thus a 


convenient basis set with which to represent the electron wave function. 


E. J. ZAK 


Zak has written a great deal on the subject of electron motion in periodic 
potentials under the influence of external electric and magnetic fields. Probably his most 
outstanding contribution to this field is the derivation and use of a Bloch-type set of 
functions which are expressible as localized Wannier functions [Ref. 14]. The 


representation of these functions is known as the kq representation [Ref. 15], and was 


derived using the fact that finite vector translations in real and reciprocal space in a lattice 
form a complete set of commuting observables. The difference between his set of 
functions and normal Bloch-type functions used to solve solid state problems is that the 
former is a function of two continuous variables k and q, while the latter are function of 
a discrete band index nand continuous wave vector k [Ref. 15]. He was an outspoken 
Critic, however, of the idea of a Stark ladder in a solid [Ref. 16]. His main criticism was 
that Wannier’s derivation was an approximation using perturbation on a one band 
Hamiltonian [Ref. 8], and, using the exact time independent Schrodinger equation, Zak 


showed the integral term of (2.19) was, in fact, an arbitrary constant, thus invalidating the 


idea of a ladder of energy levels. Wannier’s reply [Ref. 17] was to show that the integral 
term was valid, and that Stark ladders were a valid model. Zak’s reply to Wannier [Ref. 
18] was that Wannier still uses an approximate method, and his kg approach is correct. 
Finally, Zak [Ref. 19] claims, with a more robust proof, that the energy spectrum for a 
Bloch electron in an electric field is continuous. The author, thus far, has not been able to 
determine which approach is correct. Multi-band Stark levels are the electrical analog to 
Landau levels in a magnetic field, and ultimately must be proved or disproved by 


experimentalists. 


F. A. RABINOVITCH 


Rabinovitch’s first comment [Ref. 20] on the effect of an electric field on an 
electron in 2 solid concerned the translational symmetry argument in finite crystals. 
Rabinovitch [Ref. 20] showed Born-von Karman periodic boundary conditions for finite 
crystals in an electric field are incompatible with the Schrodinger equation; that is, they 
do not lead to the expected Stark ladder result for the energy levels. He [Ref. 20] also 
showed any boundary condition compatible with the Schrodinger equation breaks 
translational symmetry, and prevents the use of Bloch electron functions. His other 
works are collaborations with Zak [Refs. 21-22]. In [Ref. 22], using Zak’s kq 
representation for the Schrodinger equation, and using an identical single band 
approximation as Houston, Rabinovitch and Zak were able to derive the same result as 
Houston (2.15) for the electron wave function. They argued, since this single band 
approximation ignores other bands and interband coupling, Houston’s solution is 
unphysical and meaningless. Also, they proposed that an electron composed of the wave 


function in (2.15) would only oscillate with time faster than 2xh/eFa , and thus Bloch 


oscillations cannot exist. 
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G. J. N. CHURCHILL AND F. E. HOLMSTROM 


In a series of three papers [Refs. 23-25], Churchill and Holmstrom comment on 
the subject of Bloch oscillations and electron motion in periodic potentials subject to a 
uniform electric field. They agree with Rabinovitch and Zak’s criticism’s of this 
phenomenon [Ref. 22], and using Ehrenfest’s theorem, attempt to show why Bloch 


oscillations could not occur. They consider the following Hamiltonian of the system: 


H = p?/2m+V(x) + eFx (2.19) 


where A is a constant which does not affect the spatial distribution of the lattice potential 


V(x). They attempt to apply Ehrenfest’s theorem to (2.19) [Ref. 23] and derive the 


acceleration for the electron: 
d?{x)/dt? = {-V[AV(x)+ eFx]) = -eF +{-VV) (2.20) 


By showing a negative (non-zero) constant term (—ef’) in the acceleration, Churchill and 


Holmstrom claim [Ref. 23] the velocity must decrease with time, and thus, Bloch 
oscillations cannot occur. Churchill and Holmstrom did not show the derivation of this 
result. The author of this thesis believes this result may be incorrect, for they did not 
show that Ehrenfest’s theorem was applied to a wave function in a Bloch state (2.1). 


Also, they did not specify the periodicity of the potential V(x) =V(x +a). Their result 


appears so general as to apply to any particle accelerated by an electric field in any 
general potential, and it can be shown there are many instances where (2.20) is incorrect, 
e.g., a proton accelerated by an electric field under the influence of gravity. 

In another paper [Ref. 24], Churchill and Holmstrom build the electron wave 
function using eigenstates assuming, a priori, that the energy states are in the form of a 
Stark ladder. They show [Ref. 21] that their Bloch-type wave functions do not approach 


conventional zero field (F = 0) Bloch states. Later analysis by Krieger and Iafrate [Ref. 
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23] show this to be the fault with their analysis and later conclusions. (Note: Other 
analyses by other researchers of Bloch oscillations show that as F — OQ, the electron 
remains in a Bloch state [Refs. 8-11, 14-15, 21-22, 26]) Their final conclusion, again, is 
that Bloch oscillations do not exist, but that “... accelerated Bloch states can be looked 
upon as either a type of standing wave or produced as an interference of a stream of 
electrons moving in the +x direction and a stream moving in the —x direction.” [Ref. 24]. 
In 1991 [Ref. 25], they attempt to refute Krieger and Iafrate’s criticisms of previous work 
[Refs. 23-24], and show that Bloch oscillations exist in free space. Based on their proof 
[Ref. 25], they conclude that Bloch oscillations cannot exist in a solid if they exist in free 


space. The author of this paper believes their conclusion is incorrect. 


H. J. B. KRIEGER AND G. J. IAFRATE 


In 1985, Krieger and Iafrate [Ref. 26] summed up the criticisms [Refs. 14-16, 20- 
25] of the previous analytical solutions to Bloch oscillations by Zak, Rabinovitch, 
Churchill and Holmstrom. According to Krieger and Iafrate [Ref. 26], the criticisms can 


be summarized as follows: 


1. The energy eigenvalues of the time independent Schrodinger equation are not 
quantized but are continuous with all values of E) allowed. 


2. Since the Hamiltonian is not periodic on the boundaries of a (finite) crystal, it 
is not clear that one can employ the crystal momentum representation of Houston 
functions since Bloch functions are periodic on the boundary, i.e., a superposition 
of Bloch functions to represent the wave function w will automatically yield a 


w which is periodic on the boundary, but the solution of the time-dependent 
Schrodinger equation, including the non-periodic scalar potential, @ = -eEx may 
not have that property. 


3. The crystal momentum representation of the operator x which enters in the 
calculation may not be well defined because x,, cannot be represented as a 


linear combination of Bloch states, i-e., {|x,,|° dt diverges as the crystal 
approaches the infinite extent in the x direction. 


Krieger and Iafrate approach the solution of the time evolution of electrons in an electric 
field slightly different. Instead of using the familiar scalar potential for the electric field, 
the chose [Ref. 26] to represent it as a vector potential which is a function of time. The 


time dependent Schrodinger equation as solved by Krieger and Iafrate is: 


Hty(x,t) = eo Giey +V(x) 


L 


W(x, t) = pd (2.21) 





where A = f&(t')dt’, and &(t)is an electric field turned on at time ¢ = 0. 


The eigenfunctions [Ref. 26] are solved from following equation: 


(2. (¢ Z c)A)° . ven} Pyaneenpnney: (2.22) 


= 





Using a gauge transformation [Ref. 26], the solution of the eigenfunctions are: 


p/(x,t) = eo, , (x), (2.23) 


Expanding the wavefunction (x,t) as a sum of eigenfunctions in Eq (2.24), and solving 
for the time variation of the expansion coefficients, a, (t), of the eigenfunctions, Krieger 


and Jafrate [Ref. 26] reach the same result as Houston [Ref. 7] did: 


F(t) 


EF pen kt"))men' (RCE) 
ih 


a(t) = > Ay (LX nn (K(E))e (2.24) 


Zak [Ref. 27] criticized this formulation, and claimed this solution was incorrect due to 
the wave function repeating itself at every time on every lattice cell. He concludes the 
basis set (2.24) is incorrect, and that Krieger and Iafrate moved the problems of spatial 
periodicity into the time domain. Krieger and [Jafrate’s rebuttal [Ref. 28], agree that their 
basis states are periodic in time at lattice sites, but disagree with his conclusions 


concerning this fact. They show for any lattice site that the wave function repeats itself at 
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the Bloch repeat time, not just any time as Zak [Ref. 27] assumes. Also, they show in 
Zak’s derivation that he actuaily uses a different Hamiltonian than they do [Ref. 28], and 
seem to overcome previous objections with boundary conditions [Refs. 14-16, 20-25] 


using their expansion on basis states based on the vector potential derivation. 


i. L. ESAKI AND R. TSU 


Esaki and Tsu, in 1970, predicted that Bloch oscillations might be seen in a 
structure called a “semiconductor superlattice” [Ref. 3]. Given the state of art in 
semiconductor fabrication, it would be possible to grow layers of different material in a 
way that a periodic band structure would be created. Given the probability of Zener 
interband tunneling, Esaki and Tsu [Ref. 3] said that if the scattering times were long, 
Terahertz (THz) oscillations would be possible, and the electron would undergo 
reflections at miniband boundaries [Brillouin Zone]. The typical numbers needed, as 
estimated by Esaki and Tsu for frequency v, = 250 GHz, F =10°V/cm, and a= 10 


nm, with the scattering time greater than four picoseconds. This article was the first in 
which the possibility was raised of creating a siructure in which radiation from an 


electron undergoing Bloch oscillations might be obtained. 


J. D. EMIN AND C, ¥. HART 


In 1988, Emin and Hart [Ref. 29] commented on the time evolution of electrons in 
a periodic potential under the influence of an external electric field. Their novel approach 
showed that the applied electric field can be broken up as a sum of a periodic function 
(sawtooth), and non-periodic (ladder), and that in their multiband solution, the ladder 
forms the basis of the Stark ladder of multiband energy eigenstates [Ref. 29]. Using a 
modification of Krieger and lafrate’s approach [Ref. 26], they derive an equation of 


motion where their Bloch eigenstates are a function of the periodic portion of the electric 
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field [Ref. 29]. Since then, the Emin-Hart approach has been shown by several authors 
to be incorrect [Refs. 30-32]. 


K. M. LUBAN AND A. BOUCHARD 


In 1993, Luban and Bouchard [Ref. 5], using highly accurate numerical 
techniques, simulated the motion of an electron wavepacket in a semiconductor 
superlattice under the influence of an electric field. The advantage of their approach was 
the numerical simulation contained all powers of the Hamiltonian containing all energy 
bands, instead of single band Hamiltonian approximations used by previous analytical 
approaches. Coinciding with recent four wavemixing experiments which detected Bloch 
radiation for the first time from superlattices [Ref. 4], Luban and Bouchard’s numerical 
simulations clearly demonstrated (absent scattering events) that an electron will undergo 
Bloch oscillations in a semiconductor superlattice, and for the first time, showed that 
Bloch oscillations are a bona-fide aspect of the dynamics of independent electrons in a 


solid. 





Ill. NUMERICAL METHODS 


A. TIME DEPENDENT SCHRODINGER EQUATION (TDSE) 


Solving the Time Dependent Schrodinger Equation is essential to the modeling oz 


Bloch oscillations. The general one dimensional Schrodinger equation is: 


Hf FW, aca aL pe ney (3.1) 
2m ax? 


where w(x,7) 1s the wave function of the electron, and V(x, t)1is the potential energy 
function of the electron. We will consider a system in which V(x, 7) is independent of 


time, thus making V a function of position only. The Hamiltonian of a system is defined 
as H = T +V, where T is the kinetic energy operator, and V is the potential seen by the 


electron. The Schrodinger equation can be cast in terms of the Hamiltonian operator. If 


angie ho : 
the momentum operator p is given by Pa then the Schrodingez equation can be cast 
1 Ox 


a} 


in the form of an operator equation, with H = PV. The time dependent Schrodinger 
m 
equation reads: 
dw (x,t 
Hy(x,t) = pee (3.2) 


To solve this numerically, one must discretize (3.2). We will discuss the discretization 


process in the following sub-sections. 


1, Discretization of the Wave Function 


First, we represent the continous variable x as a discrete variable, x = x, + nAx, 
where Ax is the mesh size, x, is the starting value, and nis any integer. Also, we 
represent time as a discrete variable, t = t, + mAt , where At is the time mesh size, f, is 
the initial time, and m is any integer. The initial points, x,,7, are both set equal to zero 


to simplify the scheme. Also, this form of discretization for time allows the efficient use 
of a loop in a computer program. The form of our sampled wave function is a column 
vector, with each sample being taken at a mesh point, such that for any time ¢ , the wave 


function will look like: 


(Ax) 
tp(2Ax) 
: (3.3) 


where JN is the total number of mesh points in the discretization. Using a loop over time, 


the computer program generates a new w vector for each time step. 


2. Discretization of the Hamiltonian 


— — h°> ae 
The Hamiltonian has two terms: the kinetic energy operator T = — ae and 
m Ox 


the potential V(x). We must find a suitable discrete representation for each term. The 


discrete representation of the potential is easy. 


Being a function of x, it takes on the same form as the wave function, i.e., a column 


vector such that: 


V(Ax) 
V(2ax) 
V(x) = (3.4) 
V((N - 1)Ax)| 
V(NAx) 


The discretization of the kinetic energy operator is more involved. A suitable discrete 
representation for the second spatial derivative must be found. To approximate it, we 


look at the Taylor series expansion for arbitrary functions f(x +Ax) and f(x - Ax): 


d d 
f(x +Ax) = f(x x)+= i oO as = SIO) net. (3.5) 


and 


1 LO) ns salah entes) 


ama ie Ax” +... (3.6) 


eg) 


Ignoring terms higher than order two, and taking the sum of (3.5) and (3.6), we find that 


d*f(x)_ f(x+Ax)+ f(x-Ax)~2f(x) (3.7) 
dx? Ax? 


Letting the total derivative go to a partial derivative in x, and letting w be our arbitrary 
function f , the kinetic energy operator T acting on w for an arbitrary mesh point 


nAx is: 





Z 
h 
-——> w(nAx) = - 5 [w((n +1)Ax) + w((n -1)Axr)-2y(ndx)]. (3.8) 
m9 2mAx 
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Therefore, the Hamiltonian acting on w at an arbitrary mesh point nAx has the form: 





Hy(ndx) = - [w((n + 1)Ax) + y((2 - 1)Ax) - 2y(ndx)] + V(ndx)y(ndz) (3.9) 


2 
2mAx 


Tf our mesh runs from n = 0,1,2...N , in solving the the Hamiltonian, it arises that values 


for y(-1Ax) and w((N +1)Ax) are needed. These exist outside the mesh system we are 


looking at, and therefore are simply set to zero. This is tantamount to assuming the 
potential suddenly becomes infinite outside the mesh system. ‘This, as a practical 
consequence, forces one to construct the spatial system sufficiently large so that, over the 
time scales of interest, the wave functions do not reflect from these “hard walls.” 

The effect of the discretization is thus to replace a single differential equation with 


a set of N coupled equations. The form of the Hamiltonian isan N X N matrix: 


Nailed a 0 0 
| oi -2a+V «a 0 
= 0 oS, ws 0 (3.10) 
: 0 a -2a+V C 
0 ee 0 os —2a+V 


93 


where @ = — SMEs It is seen that the Hamiltonian is a tridiagonal matrix, with only 
ae? 





the main diagonal and the first upper and lower diagonals having non zero terms. To 
reduce the amount of space required for computer memory, the computer program in 
Appendix A uses a sparse representation for storing the Hamiltonian. Instead of an 


N X N matrix, we use an N X 3 matrix storing only the non-zero terms: 


0 -20+V a] 
a -2a0+V ai! 


a -2a¢V a 
a -2a+V @Q 


B. SOLUTIONS OF THE TIME DEPENDENT SCHRODINGER 
EQUATION 


Equation (3.2) is a first order differential equation of the time variable. The 
solution of this first order differential equation has a unique solution, given in the 
Schrodinger representation by: 


ifft 


w(x, t) =e ® y(x,0) Gal) 


To perform practical computations of w(x,t), we must know two things. First, the initial 


State of the wave function must be known or guessed for time t = 0. Secondly, an 
iss 
approximation fore ” must be made. This value is not only complex, but it is also 


Strictly unitary, and any approximation must maintain unitarity for it to be considered 
valid. Unitarity means that the normalization of the wave function is constant at all 


times, and the wave function for any time ¢ must satisfy the following integral equation: 


fy (x%,.w(a, tdx =1 (3.13) 


In the following sub-sections, we will discuss approximation methods and their relative 


strengths and weaknesses. 


1. Power Series Expansion 


_ittt 
Numerically calculating the time evolution operator e * is not easy. One could 


use a DOWer Series expansion to represent the operator. For an arbitrary operator A, the 


power series expansion of e“ is: 


2 3 
A UA 

Be (3.14) 
ens 


Based on this result, the power series expansion of the time evolution operator is: 


(3.15) 


A truncation of the expansion can be implemented easily numerically, since it is only a 
finite sum of products of terms. The restriction that any approximation maintain unitarity 
and satisfy (3.13) prevents the use of it in program. Any truncation of this infinite series 
removes the property of unitarity from it, and since a computer cannot generate an 


infinite amount of terms, any solution generated from the truncation would be invalid. 

2. Magnus Expansion 

In 1954, Wilhem Magnus derived an expansion which can be used to represent the 
time evolution operator. Using the Baker-Hausdorff theorem, Magnus derived a solution 


to the first order differential equation [Ref. 33]: 


“¥( = A(t)Y(t) (3.16) 


where Y(t) is of the form: 


Y(t) =e Y(0) (3.17) 
and 02(t) is: 
©2(t) =f A(t)dt + . fA@,fAto)dolae “g (3.18) 


where the brackets indicate commutators. This result was independently derived by 


Pechukas and Light [Ref 34]. Applying the Magnus expansion requires a very 
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sophisticated numerical integration routine, and requires a great deal of error analysis and 
monitoring of the results tc ensure numerical stability and accuracy [Ref. 35]. Because of 


its complexity, this method was not chosen. 


3. Explicit Alogrithm for Solving the TDSE 


In 1991, P. B. Visscher developed an explicit numerical algorithm to solve the 
time depedent Schrodinger equation [Ref. 36]. Starting with equation (3.1), Visscher 


rewrote it in terms of solutions of the real and imaginary portions of the wave function 


such that: 
aw (x,t 
SPE 6 Hy (x,1) (3.19) 
at 
and 
aw .(x,t 
as 2 2 —Ehy, (2,1 (3.20) 


where w, (x,t) is the real portion of the wave function, and wp ,(x,7)is the imaginary 


portion of the wave function [Ref. 36]. He claims this approach is the same as calculating 
a trajectory, and thus in the discretization of each equation, each portion is defined at 


Staggered times. For the real portion of the wave function wW_(x,f), it is defined at times 


equal to t = 0, At,2At,..., while the imaginary portion is defined at times 


t = OSAL,LSAt,25At,...[Ref. 36]. The discretization of equations (3.19) and (3.20) is: 
w(x, f+ OS5At) = wW,(x,t - OSAL) + At (x,t (3.21) 
and 


ww (x,t + OSAL) = w, {x,t — OSAL) — Athy, (x, 1) (3.22) 


To start the iteration, initial values for w,(x,0) and w,(x,05Ar) must be obtained, ana 


time must start at t = OSAt [Ref. 36]. Visscher also points out that because the real and 
imaginary portions of the wave function exist at different time steps, calculating the 
overall probability density could be problematic. A set of equations which conserve 


probability for integer and half integer units of time are: 
lp(x, 0)? =, (x,t) +, (x,t + O5At)w (x,t - 05Ar) (3.23) 


for integer time intervals, and for half integer time intervals: 


Iw(x, 0) = w, (x,t + OSAt)w, (x,t - OSAL) + W, (x,t) (3.24) 


This scheme is second order accurate, and conserves probability, satisfying equation 
(3.13). The claim of Visscher is this scheme is three times faster than the actual scheme 
used by the author of this thesis (Cayley method). This paper was discovered in the 
process of the thesis write-up, and therefore this method and its author’s claims could not 
be compared to the Cayley method. If Visscher’s claim 1s true, then it would be an 


outstanding method to use, and would be an excelient model for future students to use. 
4. Cayley Hamilton Approximation 


The approximation actually used for the time evolution operator is the Cayley- 
Hamilton approximation [Refs. 37,381. It is Hermitean, and maintains unitarity, as well 
as has good properties for keeping error low. The Cayley approximation of the time 


evolution operator is given by: 


iHAt 1: 
SX 1-3 iHAt 


e a: (3.25 
1+ LiHAt >) 
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which has error on the order of (At) . This error can be seen when compared to the 


tHAt 
power series expansion of e * . Taking the first four terms in the power series 


expansion of this time evolution operator from (3.14): 


Hist Ar (HAL) _(HAty’ 
—— tt es 


a — : 3.26 
nh one 6h? aed 
: o... 1=s1 ; 
while the expansion for the approximation ei iS: 
1+3:iHAt 
WA:  iHAt (HAN@ i(HAD? 
42 = J — —— - + +5 — (327) 


1++iHAt h me Aho 


these two expansions differ beginning only in the third term, so the error will be third 


order. The Cayley approximation to the original solution is: 


— 5, lHAt 
wW(x,t+ At) = baat 


——2h yx, t 3.28 
eae ais 


Finally, the computer will iterate the following algorithm: 


14+ iH At)w(x,t+ At) =(1-S1HAt)y(x,t) ss 
(1+ 37 VAT) ( 2 


In numerical analysis, this type of solution is called an “implicit” scheme or a Crank- 


Nicholson scheme [Ref. 38]. 


C. SOLVING THE CAYLEY APPROXIMATION 


Equation (3.29) has the form: 


M'w(x,t + At) = Mip(x,!) (3.30) 


A 


where M =] + ree , and J is the identity matrix, and H isan N XN matrix. 
zZ 


Solving for the wave function at the next time step is the ultimate goal in order to model 


the time varying motion of the electron. There are several methods of solving equation 


(3.29) available, and the methods will be discussed in the following sub-sections. 
1. Solution by Inversion 
The most direct solution of equation (3.30) for w(x,t + Af) is: 


(x,t + At)=(M") My(x,2) (3.31) 


e ° e e e ® 6 e e e e e * —] e 
but it is numerically intensive. M isa tridiagonal matrix, while the inverse, (M ) , 1s 


a full matrix, and must be solved before proceeding to solve equation (3.31). The inverse 


can be obtained by several methods (Gaussian reduction, Gauss Jordan elimination, 


etc...), but the ultimate solution of equation (3.31) takes an order of N° operations to 


solve [Ref. 38] if (M ve and Mare N X N matrices. There is a more elegant method 


which takes advantage of the tridiagonal nature of M and M_, and reduces the number 
of operations to the order of NV. This is the L-U decomposition method, and will be 


discussed in detail next. 

2. L-U Decomposition 

This method of solving a linear equation was taken from Numerical Recipes in C 
[Ref. ]. Equation (3.30) has the form A-x = b, where bis an N X 1 column vector and 


is the solution of [J - 4“ Hhp(x,t). A=[/+#H]= M , and xis what we are solving 


for: w(x,t+ At). 


26 


In order to solve this by L-U decomposition, we break A up into two matrices such that: 
A-x=(L-U)-x=L-U:x)=b. (3.32) 


We solve this by back substitution, first solving for an intermediate vector y , such that 
L-y =b, and then back solve U-x = y. This involves solving a set of N linearly 


independent equations. In decomposing a tridiagonal matrix, L takes on the form: 


oeO 0 
Cri, ues 
a ey (3.33) 


Oyewa Gyavs 8 
0... 0 Ovan Onn 


where o.,,are the coefficients of the L matrix. In order to conserve computer storage 


memeory, the program in the Appendix stored the Z matrix inan N X 2, which only 


held the non-zero elements, so that 1s: 
On 0 


L= ° (3.34) 
Oye2n-1  On-1,N-1 


| Oy iN OL wy 


a 


Likewise, the U matrix was handled in a similar fashion. The form of the U matrix when 


derived from a tridiagonal matrix is: 


Bi Bi. 9 cl 0 
0 Bx» Bos 
U=|: 5 0 (3.35) 
Byawes Bye 
40 ce 0 Baw | 


where §,. are the coefficients of the U matrix. 
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Again, to conserve space, U was actually stored in an N X 2 matrix and takes almost the 


same form as L: 


(3.36) 


By -1,N-1 B N-1.N 
0 By . 
The method of solution is straight forward. First, all the coefficients of the LZ and 
U matrices must be solved. Then equation (3.52) can be solved in two steps. The first 


step is to solve for L- y = b, which can be done as follows: 


b, 
—- 
On 
1 i-1 
Wh =a [6- Sows for i = 2,3,...,N (3.37) 


This procedure is accomplished very quickly due to the tridiagonal nature of the L 


matrix. Then to solve for x, the final answer, it is solved by back substitution as follows: 


YN 
xX, eat 
"Baw 
N 
ora 2p fori=N-1N -2....,1 (3.38) 
i jets 


Once this answer for the wave function for the new time step is obtained, it is multiplied 
by M to form anew bvector, and the process is repeated to solve for the wave function 
at the next time step. This process is repeated in a loop over the time of the simulation. 


Results of the program will be discussed in the next chapter. 


28 


IV. RESULTS AND CONCLUSIONS 


The computer code in the Appendix has been under development for nearly a 
year. The solution of the time dependent Schrodinger equation using the Cayley 
approximation sounds easy to accomplish numerically, but has been challenging to 
successfully accomplish. C language was chosen due to the local support and our 
familiarity with it. Since the electron wave function is complex, solving the time 
dependent Schrodinger equation numerically required the use of a user defined structure. 
All variables were defined of type complex, with a real and imaginary part of double 
precision. The first efforts in developing the code were writing reliable functions to add, 
subtract, multiply, divide, and conjugate complex numbers. These tasks were verified 
using simple algebraic equations that the author programmed and knew the results ahead 
of time. The functions performed as expected. 

The next step was to write a routine that would perform an L-U decomposition on 
an N X N matrix. The concepts were well laid out in Numerical Recipes in C [Ref. 38], 
but tailoring it a tridiagonal matrix was more Gifficult than expected. After several 
attempts, C code was written to perform the L-U decomposition, and results of our 
decomposition routine on a general N XN matrix were verified against an identical 
matrrix using MATLAB’s L-U decomposition function. The decomposed matrices were 
identical. Finally, the back substitution routine was written to solve the time dependent 
Schrodinger equation (3.29) using the Cayley approximation. To ensure our code was 
accurate, our first model was a one dimensional free particle, with V = Q. The wave 


packet for this model, and all subsequent simulations was a normalized Gaussian wave 


x? 
packet w(x) = Ce 2° , where Cis a normalization constant, and o controls the width of 


the packet. The free particle model behaved as expected, with the wave packet spreading 
out in both the +x and —xdirections. To ensure our numerical model was correct, 
another system simulated was an electron wave packet in a harmonic oscillator potential 


V =14mqm°*x’. The response of the particle shouid be to diffuse out, and gather itself back 


jae: 





at exact time intervals of 2/0). The model again showed the wave packet evolve as 
expected, and repeat at exact time intervals of 2n/w. 

Confident that the model to numerically simulate the time dependent Schrodinger 
equation was accurate, electron wave packets in periodic potentials with an external 
electric field applied were modeled. The simulations were ran either on a Sun 
Sparcstation 20, or a Cray C-90 super-computer at Wright Patterson AFB. For problems 
modeling Bloch oscillations, the Sparcstation typically took 12 hours to run, while the 
Super-computer would take eight hours. A variety of results were generated, and a 


representative sample is presented and discussed below: 
A. NON LOCALIZED WAVE PACKET IN GENERAL POTENTIAL 


Figure (1) is a picture of an electron clearly undergoing Bloch oscillations. 





Figure (1). Electron Undergoing Bloch Oscillations 


For all simulation output, the vertical axis represents time, with the electron wave packet 
Starting at the top of the figure working its way down in time, and the horizontal axis of 
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the figure is space. The color bar at the bottom of the figure shows the relative 
probability, with zero probability being green, and maximum probability being yellow, 
with a continuous colors scale showing probabilities in between. The horizontal scale is 
10 nm, and the amplitude of oscillation being nearly that large. For tne system in Figure 


(1), the potential function is given by V = 05|cos(2x)|, with a = 1/2 nm being the repeat 
distance. The applied electric field strength is F =10° V/cm. Given the previous 
parameters, the Bloch repeat time is t, = <t- = 264 fs . The code was run for 16000 time 


steps, which generated three oscillations. {t is clear the center of mass of the electron is 
oscillating back and forth, and thus, we would expect to radiate like a classical dipole. 
No other dynamic electron effects are seen here, so, by accident, we chose a system that 
would undergo perfect Bloch oscillations instead of executing interband transitions or 
remaining in a stationary state. Without calculating the energy bands and eigenvectors 
and composing the inital state of the electron as a superposition of eigenfunctions of a 
single band (Wannier functions), the behavior could not be exactly predicted prior to 


running the code. 


oul 





5B. LOCALIZED ELECTRON IN A GENERAL POTENTIAL 


Figure (2) is a picture of an electron undergoing a “breathing mode” of oscillation. 





Figure (2). Electron in a “Breathing Mode” 


The system is exactly the same as previous, with the only difference being that the 
electron wave packet is localized in a well. We did this by reducing 0, which controls 
the width of the wave packet. Instead of executing center of mass oscillations, the 
electron wave packet “breathed” in and out at the Bloch frequency. The electron 
preferentially breathed in the direction of lower energy. The code was run for 16000 
time steps, showing three breathing modes. This phenomenon would not radiate, 
because, as seen in the picture, the average center of mass does not appreciably move 
over time. This phenomenon was actually predicted by Zener in 1933 [Ref. 2], and his 
explanation is the electron is localized near the Brillouin zone boundary, and it only 


transitions between k = +2 and k = ~= in k space by Bragg reflection, thus executing a 


a 


“breathing” in real space. 





C. ELECTRON REFLECTIONS AT BOUNDARIES 


Figure (3) is an example of electron wave packet reflection off the boundaries of 


the system investigated: 





Figure (3). Electron Reflection of Boundaries 


In the process of discretizing the time dependent Schrodinger equation, it imposes finite 
boundaries on the system. The effect of reaching a boundary, as seen in Figure (3), is 
quite dramatic. The wave reflects perfectly off the boundary, and begins to interfere with 
the rest of the wave function executing a Bloch “breathing” oscillation. After several 
reflections, the “breathing” becomes completely distorted, as the wave function folds 
over on the right side. For future simulations, a test could be developed to determine if 
the wave function was near a boundary, and the test would stop the simulation if a 


reflection occurred. 


Sie, 





D. KRONNIG-PENNY SIMULATIONS 


The Kronnig-Penny one dimensional square wave potential is an excellent model 
to simulate the conduction band edge of a general semiconductor superlattice. The 


superlattice simulated is a GaAs/Al,,,Ga,., As superlattice, with barrier height equal to 
0.243 eV. In the following sections, we will discuss the effect of varying different 
parameters, and their effect of the motion of the electron in the lattice. 


1. Non-localized Electron in a Superiattice 


Figure (4) is a picture of an electron in a superlattice potential: 





Figure (4). Non-localized Electron in a Superlattice 


The potential function is given by V = 0.243 eV for the Al ,,Ga,,As barriers, V = 0 


eV for the GaAs barriers, and a = 10 nm being the repeat distance. The AlGaAs barriers 


were 5.1 nm thick, and the GaAs barriers were 4.9 nm thick. The applied electric field 
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strength is F =10°V/cm. Given the previous parameters, the Bloch repeat time is 
Tz = 4 = 827 fs, and the amplitude of oscillation was approximately 20 nm. The 


simulation was run for 17000 time steps, giving two oscillation periods to observe. The 
electron wave packet was approximately 250 angstroms wide. It is seen part of the 
electron is “ripped away” and undergoes Bloch oscillations, and parts of it bifurcate and 
tunnel through the barriers, remaining relatively stationary. This electron is clearly ina 
superposition of a dynamic and stationary state. The electron is undergoing Bloch 
oscillations, but it could also be transistioning energy bands or remaining in the ground 
state. Without calculating energy band transition probabilities, it is impossible to tell in 


this calculation. That could and should be built into a more sophisticated model. 


2. Effect of Varying Repeat Distance 


Figure (5) is a picture of an electron in a superlattice with half the repeat distance 


of the electron in section 1: 
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Figure (5). Electron in Superlattice (Differing Repeat Distance) 
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The repeat distance was changed from the system in section 1 from 10 nm to 5 nm. 





Keeping the electric field constant, this changes the repeat time to t, = 2 = 4135 fs. 


The simulation was run for 17000 time steps, showing four oscillation periods. The 
result is much the same as the previous two before, except the electron tunnels through 


more barriers, and less of it is available for Bloch oscillations, which are faintly present. 
3. Effect of Varying Barrier Thickness 
Figure (6) shows the effect of changing the system in section 2 by reducing the 


width of the AlGaAs barriers, and increasing the GaAs well width such that the5 nm 


repeat distance is maintained: 
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Figure (6). Electron in a Superlattice (Varying Barrier Wicth) 


The electron completely tunnels through the barriers, and it appears the electron is not 
undergoing Bloch oscillations at all. This may be because the barriers are relatively thin 


(1.5 nm), and repeat very closely to one another (5 nm). 
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4. Effect of Varying Electron Mass 


Figure (7) is a picture of the effect of changing the electron mass of the system 


described in section 3: 





Figure (7). Electron in a Superlattice (Varying Electron Mass) 


The mass of the electron in the system in section 5 was changed to 0.85 times the bare 
electron mass. The effect of making the electron “lighter” was to make Bloch oscillations 
somewhat apparent, where as they were not in the previous simulation. This represents a 
physical effect. In a semiconductor, the interaction of an electron with each material can 
be modeled by changing the electron mass to an effective mass. While the actual 
effective masses of GaAs and AlGaAs were not used, the simulation is effective in 
demonstrating the effect of the effective mass approximation in semiconductors, and 
shows that semiconductors with effectively “light” electrons are ideal materials to 


stimulate Bloch oscillations where as they might not be present on other materials. 
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E. SUMMARY 


We have demonstrated that Bloch oscillations are a genuine phenomenon of an 
electron in a periodic potential. The rich variety of dynamic phenoma of Bloch electrons 
is a great source for further research. There are a wide variety of potential applications 
for stable GHz and THz oscillators for high speed integrated circuits and Rf sensors. 
More sophisticated models of this effect can be made such as including the effective mass 
variation of the layers, calculating the energy eigenvalues and vectors to start the initial 
state in a single energy band, calculating the effect of the electron motion on the potential 
using the Poisson equation, and using that potential in the Schrodinger equation (self- 
consistent Poisson-Schrodinger equation), incorporate scattering mechanisms such as 
electron-electron scattering, electron-hole elimination, electron phonon/photon scattering, 
and a variety of other mechanisms present in the solid. Including these effects adds to the 
complexity of the model and code, and would require the immense resources of a super 


computer, but accurately modeling this effect could reveal new and useful applications 


for future devices. 





APPENDIX: CAYLEY METHOD IN ANSI C 


#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#define PI 3.14159265 


struct complex { 
double r; 
double c; 


I 

struct complex Cadd(struct complex x, struct complex y); 
struct complex Csub(struct complex x, struct complex y); 
struct complex Cmult(struct complex x, struct complex y); 
struct complex Cdiv(struct complex x, struct complex y); 


Struct complex Cea(struct complex x); 


struct complex Cconj(struct complex x); 


struct complex C1, C2, M[10002][4], Ms[10002][4], psi[10002}; 

Struct complex npsi[10002], magpsi{10002], one, eye, sum, temp; 

struct complex L[10002][3], U[10002][3], y[10002], conj[ 10002]; 

struct complex zero, temp1, V{10002], b[10002], psitj 10002], magipsi[ 10062}; 
struct complex nconj[10002], nmagpsi[10902]; 


double dt,dx, hbar,m,cl, c2, tconst, sigma, e, F; 
int 1,j,N; 
long loop; 


main() 


/* Time Dependent Schroedinger Equation - using Cayley method */ 
/* L-U Decomposition program, code based on procedure in Numerical Recipes*/ 


/* in C, pp 43-44 */ 


/* Lis alpha units */ 
/* U is beta units */ 


So 


{ 


N=10001; 

m=9.1e-31; /* Mass of electron in kg */ 
dt=0.05; /* Time step in femto seconds */ 
dx=0.01; 


hbar=0.6582; /* Plancks constant divided by 2 pi in eV fs*/ 

tconst=0.03810; /* hbar squared divided by twice mass of electron eV nm2 */ 
c2=tconst/dx/dx; /* Constant that wraps up all constants*/ 

cl=0.5*dt/hbar; /* second constant that wraps up all constants */ 
sigma=4000.0; 

eae; /* Charge of electron */ 

F=1e-2; /* Field strength in Volts/nm */ 


Cli=0.0; 
@ilkc=120" cl: 


C72 7=1.0* 2; 
CZ7c=0.0; 


zero.r=0.0; 
ZE1O C=O): 


one.r=1.0; /* One - real number */ 
one.c=0.0; 


eye.r=0.0; /* Imaginary i*/ 
eye.c=1.0; 


for (i=1; 1<=N; i++) /* Initializes M matrices to zero */ 
for G=1; j<=3; j++) 


{ 
M{ij[j].r=0.0; M[1i][}].c=0.0; 
Msf1]{}]}.r=0.0; Ms[i][3].c=0.0; 
/* End j loop */ 
} /* End i loop */ 


for (i=1; i<=N; i++) 


{ 

Vii].1=0.5 *fabs((cos(2.0*i*dx)))-(e*F*i*dx); 
V{i].c=0.0: 

; 


Jj 
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tcr (i=1; i<=N; i++) /* This initializes the M matrices in sparse format */ 
li (i==1) 


{ 
M[i}[2].r=1.0; 
Mfji][2].c=-1.0*(c1* Vfi].r+2*cl *c2); 
Ms{i][2}=Cconj(M[i]{2}); 
M[i][3].r=0.0; 
M[i}[3].c=c1 *c2; 
Msfi][3]=Cconj(MIi][3]); 
/* Endhit */ 
else if G== 
{ 
M[N][1].r=0.0; 
MI[N][1].c=c1*c2; 
Ms[N][1]=Cconj(M{N][1)); 
M[N][2].r=1.0; 
M[N][2].c=-1.0*(c1* V[i].r+2*c1*c2); 
Ms[N][2}=Cconj(MIN]|2]); 
} /Endeisenry 
else 


{ 

M[ij[1].r=0.0; 

M[i][3].r=0.0; 

Mfi][1].c=ci*c2; 

M{[i][3].c=cl*c2; 

Ms{i][1]=Cconj(M[i][1)); 

Ms[i][3]=Cconj(M[i][3)); 

M[i][2].7=1.0; 

Mf[i][2].c=-1*(c1* V[i].14+2*c1*c2); 

Ms[i][2]=Coonj(Mf[iJ[2}); 

} /* Bnd else */ 
} /* End i loop*/ 


for (i=1; i<=N; 1++) 
psi[i].r=0.0; 


psi[i].c=0.0; 
; /* End i loop */ 


J 
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for (i=1; i<=N; i++) 


psi[i].r=0.6*exp(-1*(i-5001.0)* (i-5001.0)/(2*sigma)); 


psifi].c=0.0; /* Remember - make divisor float */ 
psit[i]=Ceq(psii)); 
} /* End for i */ 


for (j=1; }<=N; j++) 


conj[j]=Cconj(psifi}); 
magipsi[j]=Cmult(conj[j],psi[}]); 
} /* End } loop */ 


/*LU Decomp procedure */ 
for G=1; i<=N; i++) 
{ 
Lfij[i].r=1.0; 
Lfi}[1].c=0.0; /* Initializes Lower main diagonal */ 
} /* End for.i */ 
U[1][1J=Ceq(Ms[1]{2]); — /* Initializes first Upper value */ 
L[1][2]=Cdiv(one,U[1][1]); _ /* Initializes first lower value */ 
Lf1][2]=Cmult(L[1][2],Ms[2][1]); 
for (j=1; }<=N; j++) 
{ 
U[j}[2]=Ceq(Ms[j][3]); /* Solves next 2 Upper values */ 


temp=Cmult(L{j][2],U[j]{2)); 
Uf[j+i][1]=Csub(Ms[j+1][2], temp); 


Lfj+1][2}=Cdiv(one,U]j+1][1]); /* Calculates next Lower value */ 
Lfj+1][2]=Cmult(L[j+1][21,Ms{j+1if1)); 
/* End j loop */ 


for (loop=9; loop<=16000; loop++) 


{ 


for (1=1; 1<=N; 1++) 
{ 
sum.r=0.0; sum.c=0.0; /* This creates "b" vector */ 
if G==1) 


temp=Cmult(M[1i][2],psit[1}); 
sum=Cadd(sum,temp); 
temp1=Cmult(M[i][3],psit[21); 
sum=Cadd(sum,temp1); 
[= end if */ 
else if (i== 


{ 
temp=Cmult(M[1i}[1],psitfN-11); 
sum=Cadd(sum,temp); 
tempi=Cmult(M[1][2],ps:t{N 1); 
sum=Cadd(sum,temp1); 
} / ABOORCISE Vi) 
EISe 
i 
for (j=1; 3<=3; J++) 
f 
temp=Cmult(M/[il[j],psit{G-1)+G-1)]); 
sum=Cadd(sum,temp); 
} /* Bd | 7/ 
(anaee|Ser-/ 


} 
b[i]=Ceq(sum); 
} /* Bandi */ 


y{1]=Cdiv(b{1},L[1]{1]); 
for (i=2;1<=N; 1++) 
sum.r=0.0; sum.c=0.0; 


temp=Cmult(L{i-1][2],y{i-1]); /* Solves intermediate matnix */ 
sum=Csub(b/[i],temp); 


temp=Cdiv(one,Lfi}[1}); 
y[i]=Cmult(temp,sum); 
} /* End i loop */ 


npsi{ N]=Cdiv(y[N],U[N][1]); 
for (i=N-1;1>=151--) 


{ 


sum.r=0.0; sum.c=0.0; 


temp=Cmult(U[i][2],npsi[i+1]); /* Solves final answer */ 
sum=Csub(y[1],temp); 


temp=Cdiv(one, U[i][1]); 


npsi{ij=Cmult(temp,sum); 
/* End i loop */ 


for (j=1; j<=N; j++) 
{ 


nconj[j]=Cconj(npsi[j]); 
nmagpsi[}]=Cmult(nconj[j],npsi[j]); 
/* End j loop */ 


if doop%100==0) 
{ 
for Q=1; j<=N; j++) 
{ 
printf("%.3f ",nmagpsi[j].r); 


} 
printf(” \n"); 


for G=1; J<=N; j++) 
i 


psit[j]=Ceq(npsifj}); 
i /* End j loop */ 


at 


} /* End loop loop */ 


retum; 
} /* End program */ 
struct complex Cadd (struct complex x, struct complex y) 


{ 


struct complex Z; 


Z.C=X.C+Y.C; 
Z 1 =XeeEy 7; 
retum Z; 


} 


struct complex Csub (struct compiex x, struct complex y) 


{ 


Struct complex Z; 


ZAC NOCH GC: 
Zhen 
return Z; 


} 


struct complex Cmult (struct complex x, struct complex y) 


f 


struct complex Z; 


Z.C=X.I"Y.C + Y.I*X.C; 
Z.1=X.I*y.r - X.C*y.C; 


Tetum Zz: 
x 


2] 


struct complex Cdiv (struct complex x, struct complex y) 


{ 


Struct complex Z; 
double d; 
4S 


d=(1.0/((y.r*y.r)+(y.c*y.c))); 
z.1=d*((x.r*y.r)+(x.c*y.c)); 
z.c=d*((x.c*y.r)-(x.r*y.c)); 


return Z; 


j 


struct complex Cea (struct complex x) 


{ 


struct complex z; 


Z Taxa 
Z Coxe: 


return Z; 


} 


Struct complex Cconj(struct complex x) 


struct complex z; 


Ziel: 
Z.c=-1 Uexe: 


return Z; 


} 
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